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RAREFIED DIFFUSION BY MODEL SAMPLING 
by Morris Perl mutter 
Lewis Research Center 


SUMMARY 

The rarefied diffusion of molecules of species a, evaporated or emitted from one 
plate, through another species p of a different molecular weight enclosed between par- 
allel walls is analyzed by the method of model sampling. The diffusing species is con- 
sidered sufficiently dilute so that only collisions with nondiffusing molecules are impor- 
tant. Maxwellian collisions are assumed. The results are given over a Knudsen range 
from free molecule to continuum flow. These results fall along one curve for all molec- 
ular mass ratios when the appropriate definition of mean free path is used. The results 
agreed with the continuum solution limiting case and elementary slip theory was found 
to give reasonable approximate results. 


INTRODUCTION 

As more rarefied environments are encountered, an understanding of diffusion proc- 
esses under noncontinuum conditions becomes of greater importance. Compared to the 
continuum diffusion, rarefied diffusion has received very little study. The present anal- 
ysis considers the rarefied diffusion of molecules a evaporated or emitted from a plane 
surface (see fig. 1). These molecules diffuse through a nondiffusing gas of species p 
which can be of a different molecular weight. The diffusing molecules are absorbed or 
condensed when incident on either surface of the enclosure. The diffusing species is con- 
sidered sufficiently dilute so that only collisions with the nondiffusing molecules are im- 
portant. The nondiffusing molecules are considered to be in Maxwellian equilibrium at 
the temperature of the plates, which is constant. Maxwellian collisions are assumed. 

The results are given for the Knudsen range from free molecule to continuum flow. 

The method of solution is model sampling or Monte Carlo. A sample history of a 
diffusing molecule is generated by randomly choosing from the appropriate distributions 



• Diffusing species a 
O Nondiffusing species p 



at points of decision in the sample history. The sample molecules are scored (i.e. , var- 
ious characteristics are recorded) as they pass scoring positions at specified locations 
across the gap. These scoring results give the macroscopic flow quantities of interest. 

Previous analyses of rarefied gas problems using model sampling were carried out 
for Couette flow with heat transfer (ref. 1) and the normal shock problem (ref. 2). The 
present results give the molecular flux and the molecular density profile across the gap. 


SYMBOLS 


b 


D 

f 

K 

K{k 2 } 

L 

l 


numerical constant, 2.6595 
molecular collision parameter 


c,c, 


\l/2 


(2kT 


m. 


a/2 


coefficient of diffusion 
probability distribution function 
intermolecular force constant, ref. 6 
factor entering in eq. (30) 
distance between surfaces 

dimensionless distance across gap or inverse Knudsen number, L/A 
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M molecular mass ratio, rn^/m^ 

m mass of molecule 

N number of samples leaving surface 0 

n number density 

n^ i** 1 component of unit vector lying in collision plane of two colliding molecules 

n slip number density 

s 

Pyp vapor pressure 

p scoring position 

Q property of molecule 

R random number chosen from uniform distribution between 0 and 1 

T absolute temperature 

V, velocity of a or p molecules 

V R relative velocity 

v, v* normalized velocities of a molecules, V/C^ or V/C^, respectively 
v' dimensionless velocity after collision 

v R dimensionless relative velocity, V R /C^ 

5 $: 

v r dimensionless radial velocity component (normalized by C a ) 

Vp normalized velocity of /3 molecules, V^/C^ 

Z coordinate across gap 

z dimensionless distance across gap, Z/L 

a diffusing species 

/3 nondiffusing species 

y angle between molecule and z-axis 

6 net sum of Q at scoring position per molecule emitted 

e molecular collision parameter 

£ dimensionless path length in z direction, X/X^ 

© collision rate 

k Boltzmann's constant 

X path length to first collision in z direction 
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mean free path based on thermal velocity, C^/O 
A mean free path length to first collision in z direction 

p molecular flux 

cp angle parameter occurring in two molecule collision 

u> dimensionless molecular collision parameter 

o> m maximum value of dimensionless collision parameters above which effects of 
collisions are negligible 

Subscripts: 

x, y indicate coordinate directions 

0 evaluated at z = 0 at emitter plate 

+, - positive, negative z direction 


ANALYSIS 

Scoring to Find Average Flow Properties 

The average quantities of interest at various positions across the gap are the molec- 
ular flux /i, which is the number of diffusing molecules crossing the gap per unit area 
per unit time, and also the molecular density n ff , which is the number of diffusing mole- 
cules per unit volume. These values are obtained at the so-called scoring positions in 
the gap (see fig. 1). If some property of the molecule is designated by Q, the average 
of Q transported across the scoring position per molecule emitted from the surface 0 
can be written as 


n(Qv z ) 

n 0, + < v z> 


= ( 6 > V z, 0 ) d6 dv z , 


0 ) f (\0 )d6 dv z, 


( 1 ) 


The 6 is the sum of Q crossing the scoring position in the positive z direction minus 
the sum of Q passing the scoring position in the negative z direction for one sample 
molecule emitted from surface 0 at velocity v z q. 
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Th e v 2 q is the component of velocity in the positive z direction of the sample 
molecule leaving surface 0 normalized by the thermal velocity C^> as given in the no- 
menclature- The joint distribution at the scoring position f(6,v z q) can be written as a 
product of a conditional distribution and a marginal distribution, namely, 
f(6 1 v o)f( v z cP" (Br. Burt Rosenbaum of Lewis aided in the formulation of eq. (1). ) 

The variance of 8 is given as 


ff2(6) = ff 62f(6 ! V z,° )f(v z,° )d6 dv z,0 - ff 5f ( 6 i v z,o) f ( v z,o) d6 dv z,o] 


(2) 


By model sampling the integral (8) in equation (1) (ref. 3, ch. 2) can be approximated as 


( 6 ) « 6 



(3) 


where N is the number of molecules sampled and 6. is the net sum of Q crossing the 
scoring position for the i m sample molecule emitted from surface 0. 

As given in the aforementioned references, it is 95 percent certain that the error, 
which is the difference between the model sampling value 6 and the true value of (5 ) , 
is less than the so-called 95 percent confidence range given by 

| 6 - ( 6)1 <M®a(6) (4) 

Vn 


where 



(5) 


Specifically, for the case Q = 1, equations (1) and (3) give 






(6) 


■{■Vi 

where 6. is the number of times the i n sample molecule passes the scoring position 
p in the positive direction minus the number of times it passes p in the negative direc- 
tion. As shown in equation (6), the average gives the diffusing molecular flux passing 


5 


the scoring position divided by the diffusing molecular flux leaving the surface 0. The 
95 percent confidence range is given by equations (4) and (5). 

For the case Q = 1/v . equations (1) and (3) yield 


0,+( v z ) 0, + 


N 

= IV6. 

N Zf i,n 


( 7 ) 


where S i n is the sum of the absolute values of the reciprocals of the velocities (a sam- 
ple molecule can cross a given scoring position any number of times) with which the l 
sample molecule passes the scoring position. Inasmuch as the average molecular veloc- 
ity leaving the surface 0 is assumed to be given by 

<v z > 0 = mr 1/2 (8) 

where M is the molecular mass ratio m^/m^; and since the flux is merely + ( v z )q> 
equation (7) becomes 


n 


n 


0, + 


1 

(ttM) 1//2 



N i=l 


6. 

i,n 


By using equation (5), the variance can be seen to be 


cr 





1 

(ttM) 


1_ 

N 



which can be used to calculate the 95 percent confidence range (eq. (4)). 


( 9 ) 


( 10 ) 


Relation Between Mass Flux Leaving Surface and Vapor Pressure 

The molecular flux evaporating from the surface can be found by relating the molec- 
ular density ng + of the molecules leaving the surface to their vapor pressure Pyp. 

The number density of the molecules leaving the surface n n can be related to the 

u, + 

vapor pressure Pyp of the diffusing a molecules by the ideal gas law, 
n 0 + = •^v'p/2 /f 7’ (as given in ref. 4, Sec. 1. 6). 
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Choosing Velocity Components for Emitted Sample Molecules 


O j|c 

The molecular flux leaving the surface 0 in velocity increment d°v is given by 


2n 0, + v z / *2 N ,3 * 
— ^ exp(-v )d v 

.3/2 


(ID 


* 

where v = V/C . If expression (11) is divided by the total molecular flux leaving the 
surface ng + (v*) = ng + j\I ~n, the result is the velocity probability distribution function 


for the sample molecules leaving the surface 


f 0 (v ) =• 


2v 


- exp(-v* 2 ) 


( 12 ) 


This can be written in cylindrical coordinates as the product of 

f ( V z) =2v z ex P(~ v z 2 ) (13a) 

and 

f ( v r) = 2v r ex p(' v r 2 ) ( 13b ) 

Following the method of choosing randomly from a distribution as given in refer- 
ence 3 (ch. 7. Sec. 4) yields 


R = Z 2v* exp(-v* 2 )dv* (14) 

where R is a random number generated by the computer from a uniform distribution in 
the interval [0, 1]. This equation can be integrated to give 

v* = v,, Vm = (-In R) 1//2 (15) 

Z Zt 

so that by having the computer generate a random number which is substituted into equa- 
tion (15) we can solve for v„, the velocity component for the sample leaving the surface 
0 normalized by C^. A similar equation applies for the radial velocity component v p 
for the sample molecule leaving the surface. 
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PATH LENGTH TO COLLISION 


After the sample molecule leaves the surface, it is necessary to calculate its path 
length to first collision. The probability distribution function for path length to collision 
can be shown to be (ref. 1, p. 460). 


m = 



(16) 


where A is the distance in the z direction that the sample molecule will travel before 
collision. The mean free path in the z direction X is given by 


v z C 6 

A = -=-£ = vA a 
M q z 0 


(17) 


where A^ = C^/0. 

Now an expression for the mean collision frequency 0 has to be determined. The 
collisions between an a and a f3 molecule are assumed, for simplicity, to be Max- 
wellian in type; that is, the intermolecular force law between colliding molecules can be 
described by 


F - 



(18) 


where r is the distance between colliding molecules and K is a constant. It is shown, 
for example (ref. 5, p. 170), that the angle through which a molecule is deflected in a 
Maxwellian collision depends on the nondimensional quantity 


on 


' M 
VM + 1 


r2 m j3 


A/4 


'R 


K 


(19a) 


where large values of on correspond to small angles of deflection. Here is the rel- 
ative velocity before collision and b is the minimum separation distance between the two 
colliding molecules if they were not deflected from straight line paths by the collision. 

If the angle of deflection effected by a collision is taken as a measure of the degree of 
perturbation in the path of the diffusing molecule, then only those collisions in which the 
angles of deflection are larger than some small value need be considered significant. 
Equivalently, a maximum value of on denoted as on m may be established where colli- 
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sions with values of u> larger than o> m may be ignored. Thus, equation (19a) with o> 
replaced by a> m constitutes the definition of b max 


b 


max 



K y« 

v i m J 


(19b) 


where b max is the largest value of b for which a significant collision takes place. 

Proceeding as in reference 6 (p. 154) reveals that the mean collision frequency 0 
is given by 


0 = n, 


J 


W v 




max 


b db 



(20a) 


which can be easily integrated to yield the expression 





(20b) 


Acceptable values of a> m were found by carrying out the analysis for different values of 
until an increase in u> m did not significantly change the results. This aspect of the 
problem is discussed below and in reference 6 (p. 122). 

Substituting equation (20b) into equation (17) gives 



(21) 


If £ is defined as X/X^, equation (16) becomes 



To randomly choose a normalized distance £ to collision, proceed as before 

C = -v z In R (23) 
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The dimensionless distance or number of path lengths across the gap is then given as 


l ft ~ L//X j3 


where l ^ is the inverse of the usual definition of Knudsen number. 


( 24 ) 


SAMPLE MOLECULE VELOCITY AFTER COLLISION 


To find the velocity components of the sample molecule after collision, first pick 
velocity components for a collision partner. For Maxwellian collision the probability 

O 

that a collision partner will have a velocity v^ in d v^ is given by (ref. 1, p. 469) 


d0 

0 


= { /3 d3 V 


.3/2 


- exp 


W) 


ld 3 v. 


(25) 


where f^ is the velocity probability distribution function of the nondiffusing ft mole- 
cules. The function f^ can be written in cylindrical coordinates as 

f( V = 2v /3r ex p(“ v jgr) 0 ~ V ~ ~ (26a) 

0 ^ 9 < 2n (26b) 

S oo (26c) 


m =• 


2tt 




r. 


The target molecule velocity components can be picked using 

V/?r = (-In R r ) 1/2 (27a) 

v f3z = ( " ln R z )1/2 c0S ( R 2 ?r) (27b) 

6 = 2nR d (27c) 

Now having the sample molecule and target molecule components of velocity before colli- 
sion allows the velocity after collision as shown in references 6 (p. 163) and 5 (p. 57) to 
be calculated. Calling the velocity after collision v' enables the following to be written: 
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( 28 ) 


v i - v i = -(“) ( v R"i cos V sin 9 + v Ri CO s 2 <p) 
where i = 1, 2, and 3 correspond to the x, y, and z components, respectively, and 

(29a) 




R 


v Rx v Ry sin e - v R v Rz cos e 

y = 1 


V R ( v Ry + v Rz) 


1/2 


(29b) 


and the v^ is given by 


where 


n z = 


v Rx v Rz sin e + v R v Ry cos 6 


V- 


7~2 2 W 2 

R v Ry + v Rz / 


2 2 2 
Vo = l Vo,. + V Dtt + V T 


-(r 


R " VRx T Ry T Rz^ 


v Rx = v x ~ v jSx 


and so forth. The value of cp as shown in reference 6 (p. 158) is given by 


( 2 Y 1/4 fl 

/ 2 Y 1/2 

1 


1 - 1 + — 


v ^ l 2 

L v j) \ 

j 


where 


9 z* 7 r/2 o 2 .”1/2 

K{k 2 } = / (1 - k z sin z 0) d0 

^0 


(29c) 


(29d) 


(29e) 


(30a) 


(30b) 


The values of e and w are obtained by random selection as explained in refer- 
ence 6 (p. 161) from the relations 


e = 2nR, 


(31a) 
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(31b) 


a> = co m R 


1/2 

u> 


Although the value of u> m is theoretically infinite, as explained previously, collisions 
with values larger than some value of <x> m would have a negligible effect on the results. 
In references 2 and 6, where M = 1, a value of w m =1.5 gave good results. This 
value was used in the present analysis for M = 1. However, for M - 0. 1 and M = 0. 01, 
it was found that the larger value for u> m of 2 was needed; whereas, for M = 10, a 
value of w m =1.5 was used although a somewhat smaller value would have been ac- 
ceptable . 


MAXWELL'S CONTINUUM SOLUTION 


The continuum diffusion equation from reference 7 (p. 247) can be written as 


n 


jj. _ -Dw 


1/2 


0 ,+/ 




0 ,+ 


LC 


a 


dz 


(32a) 


where z equals Z/L. The coefficient of diffusion D is given in reference 7 (p. 247) as 


D 


a 


2A l n /3 


(M + l)Mm^ 


1 1/2 


K 


(32b) 


where is given in reference 7 (p. 239) as 2. 6595. This can be written in terms of 
the mean free path as 




: Vm> 


{A-, 


A. 


(33) 


where 


X M " 


/ m + l VAnV 

^2Vm A 1 - 5 / 


X /3 


(34) 


12 


Substituting equation (33) into equation (32) and letting l M = L/a m give for the contin- 
uum solution 






4.711 \ n 0,+/ 
l M dz 


(35) 


For the continuum solution, there is a no slip at the boundary so that, at z = 0, 
n / n 0 + = 2 and, at z = 1, n/n^ + = 0. Hence, 9(n/ng + ) / /sz = -2 and the continuum so- 
lution from equation (35) is equal to 


;i \ _ 9. 422 

l M 


(36a) 


This continuum solution is plotted in figure 2. It is evident that it is only a function of the 
number of free paths across the gap when based on a free path length (solid line). 
These results apply regardless of molecular mass ratios M and temperature of the sys- 
tem T. 


Molecular 
mass ratio, 

M = m a /mp 

Y 0.01 

A .1 
O 1 

D 10 

For all M or M = 1 (based on l M ) 

For given value of M (based on lj) 


Monte Carlo Conti n u u m 

results for solution for 

given M given M 



or 

Dimensionless gap distance, l^ = [(M+ 11/2]^!^ forgiven M 

Figure 2. - Molecular flux across gap divided by molecular 
flux from evaporating surface for various mass ratios. 


13 


The effects of varying the molecular masses on the diffusion rate can be shown by 
defining a dimensionless gap distance which is independent of the molecular masses. 
Using equations (21) and (34) results in 

/ 2 / K 

h=(“) L < S6b > 

Assuming that K, the term in the intermolecular force law (eq. (18)), is independent of 
M means equation (36a) can be written as 



9,422 ( M + 1 \ 1//2 

V 2 ) 


(36c) 


This result is plotted in figure 2 as the dashed lines for different values of M. The re- 
sults indicate that, for larger values of the mass of the diffusing molecules compared 
to the mass of the nondiffusing molecules (i.e. , increasing M), the diffusion rate across 
the channel increases, again assuming that the intermolecular force constant K is un- 
affected by the change in mass ratio. 


Elementary Slip Theory 

For conditions close to continuum, a common method of analysis is the slip approxi- 
mation. The usual continuum equations are assumed to apply in the gap but a discontinu- 
ity in density is assumed to occur at the wall; that is, n -f- n where n = 2n n and n„ 
is the density in the gas next to the wall. 

Consider some plane z located somewhere in the gas. The molecular flux of dif- 

ir 

fusing a molecules crossing the plane in the positive z direction can be written as 
(ref. 5, p. 307) 


( n < v z» + p 



(37) 


where the limits of integration are v x and v^ from -°° to +°° 
Expanding the molecular density in a Taylor series around point 
erence 5 (p. 307) results in 


and v„ from 0 to +°°. 
Zp as described in ref- 
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n(z) = n(z ) - X cos y(— \ 

P S \3zA 


where X g is some average mean free path and y is the angle between the z axis and 
the incident molecule. Then equation (37) can be integrated to obtain 


(n<V » = ^ (— + (vj) - — (—) (— ) 

z +P 2 \Vn z ) 3 \3z JAV^J 


Similarly, the molecular flux in the negative z direction is given by 


(n(V» =^(— -<V >) 

z -P 2 VVtt z / 3 y/jr\dzL 


If it is assumed that, close to continuum conditions, the mean diffusion velocity is small 
compared to the mean thermal velocity, (V z ) /(C a /Vn) « 1. The net flux at point p is 
then given by the difference between equations (39) and (40) as 


M = - 


2X„ C - 
s a on 

3 sfn dz 


which can be written in dimensionless form as 


H _ 2 ( X s\ V n 0,+ 


3 \L 


To calculate the slip, consider a cross section z adjacent to the surface at z = 0. 
In that case the mass flux coming off the wall in the positive z direction is given by 
equation (8). Then 


„ - 


n +0 C Q! 


+o v* 


The mass flux in the negative z direction as in equation (40) is now given as 


(n(V » + 

z -0 2 v? 3 v? 3z 


(43b) 
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where n g is the density of the gas adjacent to the wall. Then the net flux is given by the 
difference between equations (43a) and (43b) which yields 


<"< V Z »0 


n +0 C a n s C a x s C a 3n 
Vn 2 y/n 3 y/w dz 


By using equation (41) this can be rewritten as 

C. 


I 1 =■ 


a 


2 y/n 


(2n 0, + " n s> + f 


or 


1 ii 




(44) 


(45a) 


(45b) 


If it is assumed that the continuum solution as given by equation (35) applies in the gas 
but that the density adjacent to the wall is now n s /n Q + instead of 2, the continuum so- 
lution with slip would give from equation (35) 


M 


M 


0 ,+ 



(46) 


It can be shown that the results are antisymmetric around the center of the gap z = 1/2 
so the results will only be carried out for half the gap. Combining equations (45) and (46) 
results in the following for the molecular flux with slip: 




Jm\ 1 

9. 422 ) 


(47) 


This result is shown in figure 2 as the dashed double dot line. Similarly from equa- 
tions (45) and (46) the number density slip can be obtained as follows: 
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Molecular mass ratio, 
M = m^mp 

□ 10 

O 1 

A .1 

V .01 


Continuum solutions based 
on calculated flux 



Figure 3. - Number density of diffusing molecules across 
gap divided by number density of molecules leaving from 
evaporating surface. 


n 


n 


0 , + 


2l M -f 9. 422 
*M + 9 - 422 


(48) 


This result is shown in figure 3. 


RESULTS 

The resulting molecular flux for different values of the molecular mass ratio 
M = m Q /nip are plotted in figure 2 as a function of Z^ which is the number of mean 
free paths across the gap based on the mean free path These results were obtained 

for given values of Z M and M for acceptable values of The A^ was calculated 

from these values by equation (34). This gave the normalized distance across the chan- 
nel Ip which was used as input in the calculation (eq. (24)). When the results are 
plotted as a function of Z^, the molecular flux in figure 2 falls along one curve for all 
values of the molecular mass ratio M. This curve is shown as the solid line drawn 


17 


through the results. Also shown is Maxwell's continuum solution as a function of 
This also plots along one curve for all values of M. The rarefied solution approaches 
the continuum solution for large values of 

The effect of varying the molecular mass ratio M of the diffusing to the nondiffus- 
ing molecules on the molecular flux can be seen by keeping the dimensionless distance l 
across the gap independent of the mass ratio M. From equation (36b) it is evident that 
if the results are replotted as a function of l p which is the dimensionless distance 
across the gap independent of the molecular mass ratio, the effect of the molecular mass 
ratio on the molecular flux transport across the gap can be observed. These results, 
which are shown as dashed lines, indicate much higher diffusion rates for heavier mole- 
cules. In a similar manner the effect of the temperature of the system can be found, 
since the number of mean free paths across the gap can be related to the temperature of 
the system by (Tg/Tj)? - ^ = Z^/Zg (eq. (36b)). Thus, the square root of a temperature 
ratio increase equals the decrease in the path length ratio across the gap, thereby giving 
an increasing molecular flux as the system temperature increases. Also shown in fig- 
ure 2 is the simple slip solution as a function of 1-^. The simple relation given by 
equation (47) gives a good approximation to the reported numerical result. 

In figure 3 the number density of the diffusing a species at various positions 
across the gap are shown. They are plotted for given values of Z The density pro- 
files are equivalent for the same values of Z^ for all molecular mass ratios. This is 
similar to the result found for the molecular flux. The curves were found to be odd 
functions about the center of the gap. Plotted as solid lines in figure 3 are the linear 
density results based on the calculated Monte Carlo molecular flux given in figure 2 
using the continuum solution of equation (46). Shown as dashed lines in figure 3 are the 
slip solutions given by equation (48). The slip solution gives good results at the higher 
values of Z-^, but the agreement is poorer for smaller values of Z^. 

The amount of time needed to run a given number of sample histories depended 
strongly on the number of collisions for each sample, so that as the number of mean 
free paths across the gap increased the computing time increased. Using an IBM 7094 
for M = 0. 1 and 10 000 samples, the results for Z^ =1.0 took 2 minutes; for 
^M = 10, 10 minutes; and for Z^ - 30, 25 minutes. The parameters that shorten the 
mean free path will accordingly lengthen the running time. As can be seen from equa- 
tion (34), both M and to -will have an effect on the value of Z-^. 

The 95 percent confidence range as given by equation (4) was also calculated. The 
value of a {jJ-/ iIq + ) was about 0. 5 so that for 10 000 samples the 95 percent confidence 
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range for m/Mq + is about 0. 01. For the number density results the variance was 

higher and about 5. This gave a 95 percent confidence interval for n/n n of about 0. 1. 

u, + 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, February 13, 1970, 

129-01. 
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